1 DATA

#############################
#PERFORMANCE 
#############################

data_PERF <- import_data(dataset = "DATACOMPLET_PERF.csv", 
                         trait = "performance",
                         remove_testenvt = c("Grape","GF"), 
                         remove_pop = c("WT3"), 
                         remove_rate = NA)
## [1] "Data (541 and 602 tubes for the first and third generation respectively) where i) the number of eggs was NA (5 and 0 tubes for the first and third generation respectively); or ii) the number of adults  was NA (0 and 0 tubes for the first and third generation respectively); or iii) the number of eggs was zero -Emergence rate = NaN- (99 and 0 tubes for the first and third generation respectively); or iv) the number of adults was higher than the initial number of eggs (50 and 1 tubes for the first and third generation respectively) were not removed."
head(data_PERF)
##    Generation Experiment Original_environment   Population    Date_P Date_C_O
## 8          G0 Plasticity           Blackberry Blackberry45 3/10/2018  4/10/18
## 13         G0 Plasticity           Blackberry Blackberry45 3/10/2018  4/10/18
## 14         G0 Plasticity           Blackberry Blackberry45 3/10/2018  4/10/18
## 16         G0 Plasticity           Blackberry Blackberry45 3/10/2018  4/10/18
## 20         G0 Plasticity           Blackberry Blackberry45 3/10/2018  4/10/18
## 24         G0 Plasticity           Blackberry Blackberry45 3/10/2018  4/10/18
##      Date_C_A Row Column Rack Test_environment Nb_eggs Obs_O Nb_adults Obs_A
## 8  19/10/2018  L1     C3    1           Cherry       3    LO         2    LO
## 13 19/10/2018  L1     C8    1       Blackberry       1    LO         2    LO
## 14 19/10/2018  L1     C9    1           Cherry       2    LO         1    LO
## 16 19/10/2018  L2     C1    1       Blackberry       1    LO         0    LO
## 20 19/10/2018  L2     C5    1       Strawberry       2    LO         2    LO
## 24 19/10/2018  L2     C9    1           Cherry       1    LO         1    LO
##    EggScore EggScoreFive EggScoreSmall SA IndicG0 IndicG2 SAIndicG0
## 8         1            1             1  0       1       0         0
## 13        1            1             1  1       1       0         1
## 14        1            1             1  0       1       0         0
## 16        1            1             1  1       1       0         1
## 20        1            1             1  0       1       0         0
## 24        1            1             1  0       1       0         0
##                fruit_hab             fruit_hab_ng     fruit_gen       hab_gen
## 8      Blackberry_Cherry     Blackberry_Cherry_G0 Blackberry_G0     Cherry_G0
## 13 Blackberry_Blackberry Blackberry_Blackberry_G0 Blackberry_G0 Blackberry_G0
## 14     Blackberry_Cherry     Blackberry_Cherry_G0 Blackberry_G0     Cherry_G0
## 16 Blackberry_Blackberry Blackberry_Blackberry_G0 Blackberry_G0 Blackberry_G0
## 20 Blackberry_Strawberry Blackberry_Strawberry_G0 Blackberry_G0 Strawberry_G0
## 24     Blackberry_Cherry     Blackberry_Cherry_G0 Blackberry_G0     Cherry_G0
##            pop_gen      Rate
## 8  Blackberry45_G0 0.6666667
## 13 Blackberry45_G0 2.0000000
## 14 Blackberry45_G0 0.5000000
## 16 Blackberry45_G0 0.0000000
## 20 Blackberry45_G0 1.0000000
## 24 Blackberry45_G0 1.0000000
tapply(data_PERF$Nb_adults,list(data_PERF$Original_environment,data_PERF$Generation),length)
##             G0  G2
## Blackberry 169 320
## Cherry     233 143
## Strawberry  35 139
#############################
#EMERGENCE RATE
#############################

data_PERF_Rate <- import_data(dataset = "DATACOMPLET_PERF.csv", 
                         trait = "performance",
                         remove_testenvt = c("Grape","GF"), 
                         remove_pop = c("WT3"), 
                         remove_rate = TRUE)
## [1] "Data (541 and 602 tubes for the first and third generation respectively) where i) the number of eggs was NA (5 and 0 tubes for the first and third generation respectively); or ii) the number of adults  was NA (0 and 0 tubes for the first and third generation respectively); or iii) the number of eggs was zero -Emergence rate = NaN- (99 and 0 tubes for the first and third generation respectively); or iv) the number of adults was higher than the number of eggs (50 and 1 tubes for the first and third generation respectively) were removed."
head(data_PERF_Rate)
##    Generation Experiment Original_environment   Population    Date_P Date_C_O
## 8          G0 Plasticity           Blackberry Blackberry45 3/10/2018  4/10/18
## 14         G0 Plasticity           Blackberry Blackberry45 3/10/2018  4/10/18
## 16         G0 Plasticity           Blackberry Blackberry45 3/10/2018  4/10/18
## 20         G0 Plasticity           Blackberry Blackberry45 3/10/2018  4/10/18
## 24         G0 Plasticity           Blackberry Blackberry45 3/10/2018  4/10/18
## 33         G0 Plasticity           Blackberry Blackberry45 3/10/2018  4/10/18
##      Date_C_A Row Column Rack Test_environment Nb_eggs Obs_O Nb_adults Obs_A
## 8  19/10/2018  L1     C3    1           Cherry       3    LO         2    LO
## 14 19/10/2018  L1     C9    1           Cherry       2    LO         1    LO
## 16 19/10/2018  L2     C1    1       Blackberry       1    LO         0    LO
## 20 19/10/2018  L2     C5    1       Strawberry       2    LO         2    LO
## 24 19/10/2018  L2     C9    1           Cherry       1    LO         1    LO
## 33 19/10/2018  L3     C8    1       Blackberry       2    LO         2    LO
##    EggScore EggScoreFive EggScoreSmall SA IndicG0 IndicG2 SAIndicG0
## 8         1            1             1  0       1       0         0
## 14        1            1             1  0       1       0         0
## 16        1            1             1  1       1       0         1
## 20        1            1             1  0       1       0         0
## 24        1            1             1  0       1       0         0
## 33        1            1             1  1       1       0         1
##                fruit_hab             fruit_hab_ng     fruit_gen       hab_gen
## 8      Blackberry_Cherry     Blackberry_Cherry_G0 Blackberry_G0     Cherry_G0
## 14     Blackberry_Cherry     Blackberry_Cherry_G0 Blackberry_G0     Cherry_G0
## 16 Blackberry_Blackberry Blackberry_Blackberry_G0 Blackberry_G0 Blackberry_G0
## 20 Blackberry_Strawberry Blackberry_Strawberry_G0 Blackberry_G0 Strawberry_G0
## 24     Blackberry_Cherry     Blackberry_Cherry_G0 Blackberry_G0     Cherry_G0
## 33 Blackberry_Blackberry Blackberry_Blackberry_G0 Blackberry_G0 Blackberry_G0
##            pop_gen      Rate
## 8  Blackberry45_G0 0.6666667
## 14 Blackberry45_G0 0.5000000
## 16 Blackberry45_G0 0.0000000
## 20 Blackberry45_G0 1.0000000
## 24 Blackberry45_G0 1.0000000
## 33 Blackberry45_G0 1.0000000
tapply(data_PERF_Rate$Nb_adults,list(data_PERF_Rate$Original_environment,
                                     data_PERF_Rate$Generation),length)
##             G0  G2
## Blackberry 129 319
## Cherry     224 143
## Strawberry  34 139
tapply(data_PERF_Rate$Rate,data_PERF_Rate$Generation,mean)
##        G0        G2 
## 0.3593365 0.1857181
## To test Rate higher than 1
# data_PERF_Rate <- import_data(dataset = "DATACOMPLET_PERF.csv", 
#                          trait = "performance",
#                          remove_testenvt = c("Grape","GF"), 
#                          remove_pop = c("WT3"), 
#                          remove_rate = NA)
# #Replace Rate>1 by 1
# data_PERF_Rate$Rate[data_PERF_Rate$Nb_adults>=data_PERF_Rate$Nb_eggs] <- 1

###########################
#PREFERENCE 
###########################
data_PREF <- import_data(dataset = "DATACOMPLET_PREF.csv", 
                         trait = "preference",
                         remove_testenvt = NA, 
                         remove_pop = c("WT3"), 
                         remove_rate = NA)

head(data_PREF)
##   Generation Experiment BoxID     Date_P Original_environment   Population Line
## 1         G0 Plasticity   860 19/09/2018           Blackberry Blackberry33    1
## 2         G0 Plasticity   860 19/09/2018           Blackberry Blackberry33    1
## 3         G0 Plasticity   860 19/09/2018           Blackberry Blackberry33    1
## 4         G0 Plasticity   860 19/09/2018           Blackberry Blackberry33    1
## 5         G0 Plasticity   860 19/09/2018           Blackberry Blackberry33    2
## 6         G0 Plasticity   860 19/09/2018           Blackberry Blackberry33    2
##   Column Test_environment Nb_eggs   Date_C_O Obs_O SA IndicG0 IndicG2 SAIndicG0
## 1      1        Cranberry       0 13/12/2018    CD  0       1       0         0
## 2      2              Fig       0 13/12/2018    CD  0       1       0         0
## 3      3        Raspberry       0 13/12/2018    CD  0       1       0         0
## 4      4         Rosehips       0 13/12/2018    CD  0       1       0         0
## 5      1             Kiwi       0 13/12/2018    CD  0       1       0         0
## 6      2       Strawberry       1 13/12/2018    CD  0       1       0         0
##               fruit_hab             fruit_hab_ng     fruit_gen       hab_gen
## 1  Blackberry_Cranberry  Blackberry_Cranberry_G0 Blackberry_G0  Cranberry_G0
## 2        Blackberry_Fig        Blackberry_Fig_G0 Blackberry_G0        Fig_G0
## 3  Blackberry_Raspberry  Blackberry_Raspberry_G0 Blackberry_G0  Raspberry_G0
## 4   Blackberry_Rosehips   Blackberry_Rosehips_G0 Blackberry_G0   Rosehips_G0
## 5       Blackberry_Kiwi       Blackberry_Kiwi_G0 Blackberry_G0       Kiwi_G0
## 6 Blackberry_Strawberry Blackberry_Strawberry_G0 Blackberry_G0 Strawberry_G0
##           pop_gen
## 1 Blackberry33_G0
## 2 Blackberry33_G0
## 3 Blackberry33_G0
## 4 Blackberry33_G0
## 5 Blackberry33_G0
## 6 Blackberry33_G0
tapply(data_PREF$Nb_eggs,list(data_PREF$Original_environment,data_PREF$Generation),length)
##              G0   G2
## Blackberry  696 1176
## Cherry     1200  624
## Strawberry  252  480
###########################
#PREFERENCE 3 fruits
###########################
levels_test<-levels(data_PREF$Test_environment)
levels_original<-levels(data_PREF$Original_environment)


data_PREF_three <- import_data(dataset = "DATACOMPLET_PREF.csv", 
                         trait = "preference",
                         remove_testenvt = usefun::outersect(levels_test, 
                                                             levels_original), 
                         remove_pop = c("WT3"), 
                         remove_rate = NA)

head(data_PREF_three)
##    Generation Experiment BoxID     Date_P Original_environment   Population
## 6          G0 Plasticity   860 19/09/2018           Blackberry Blackberry33
## 9          G0 Plasticity   860 19/09/2018           Blackberry Blackberry33
## 11         G0 Plasticity   860 19/09/2018           Blackberry Blackberry33
## 15         G0 Plasticity   884 19/09/2018           Blackberry Blackberry33
## 16         G0 Plasticity   884 19/09/2018           Blackberry Blackberry33
## 20         G0 Plasticity   884 19/09/2018           Blackberry Blackberry33
##    Line Column Test_environment Nb_eggs   Date_C_O Obs_O SA IndicG0 IndicG2
## 6     2      2       Strawberry       1 13/12/2018    CD  0       1       0
## 9     3      1       Blackberry       0 13/12/2018    CD  1       1       0
## 11    3      3           Cherry       0 13/12/2018    CD  0       1       0
## 15    1      3           Cherry       0 13/12/2018    CD  0       1       0
## 16    1      4       Strawberry       0 13/12/2018    CD  0       1       0
## 20    2      4       Blackberry       0 13/12/2018    CD  1       1       0
##    SAIndicG0             fruit_hab             fruit_hab_ng     fruit_gen
## 6          0 Blackberry_Strawberry Blackberry_Strawberry_G0 Blackberry_G0
## 9          1 Blackberry_Blackberry Blackberry_Blackberry_G0 Blackberry_G0
## 11         0     Blackberry_Cherry     Blackberry_Cherry_G0 Blackberry_G0
## 15         0     Blackberry_Cherry     Blackberry_Cherry_G0 Blackberry_G0
## 16         0 Blackberry_Strawberry Blackberry_Strawberry_G0 Blackberry_G0
## 20         1 Blackberry_Blackberry Blackberry_Blackberry_G0 Blackberry_G0
##          hab_gen         pop_gen
## 6  Strawberry_G0 Blackberry33_G0
## 9  Blackberry_G0 Blackberry33_G0
## 11     Cherry_G0 Blackberry33_G0
## 15     Cherry_G0 Blackberry33_G0
## 16 Strawberry_G0 Blackberry33_G0
## 20 Blackberry_G0 Blackberry33_G0
tapply(data_PREF_three$Nb_eggs,list(data_PREF_three$Original_environment,
                                    data_PREF_three$Generation),length)
##             G0  G2
## Blackberry 174 294
## Cherry     300 156
## Strawberry  63 120
resume_design<-tapply(as.factor(data_PREF_three$BoxID),list(data_PREF_three$Population,
                                  data_PREF_three$Generation),length)
mean(resume_design[,1], na.rm = TRUE)/3
## [1] 7.782609
mean(resume_design[,2], na.rm = TRUE)/3
## [1] 7.916667
resume_design<-tapply(data_PERF$Nb_eggs,list(data_PERF$Population,
                                  data_PERF$Generation),length)
mean(resume_design[,1], na.rm = TRUE)/3
## [1] 6.069444
mean(resume_design[,2], na.rm = TRUE)/3
## [1] 8.026667
dim(data_PREF_three)
## [1] 1107   21
dim(data_PREF)
## [1] 4428   21
head(data_PREF)
##   Generation Experiment BoxID     Date_P Original_environment   Population Line
## 1         G0 Plasticity   860 19/09/2018           Blackberry Blackberry33    1
## 2         G0 Plasticity   860 19/09/2018           Blackberry Blackberry33    1
## 3         G0 Plasticity   860 19/09/2018           Blackberry Blackberry33    1
## 4         G0 Plasticity   860 19/09/2018           Blackberry Blackberry33    1
## 5         G0 Plasticity   860 19/09/2018           Blackberry Blackberry33    2
## 6         G0 Plasticity   860 19/09/2018           Blackberry Blackberry33    2
##   Column Test_environment Nb_eggs   Date_C_O Obs_O SA IndicG0 IndicG2 SAIndicG0
## 1      1        Cranberry       0 13/12/2018    CD  0       1       0         0
## 2      2              Fig       0 13/12/2018    CD  0       1       0         0
## 3      3        Raspberry       0 13/12/2018    CD  0       1       0         0
## 4      4         Rosehips       0 13/12/2018    CD  0       1       0         0
## 5      1             Kiwi       0 13/12/2018    CD  0       1       0         0
## 6      2       Strawberry       1 13/12/2018    CD  0       1       0         0
##               fruit_hab             fruit_hab_ng     fruit_gen       hab_gen
## 1  Blackberry_Cranberry  Blackberry_Cranberry_G0 Blackberry_G0  Cranberry_G0
## 2        Blackberry_Fig        Blackberry_Fig_G0 Blackberry_G0        Fig_G0
## 3  Blackberry_Raspberry  Blackberry_Raspberry_G0 Blackberry_G0  Raspberry_G0
## 4   Blackberry_Rosehips   Blackberry_Rosehips_G0 Blackberry_G0   Rosehips_G0
## 5       Blackberry_Kiwi       Blackberry_Kiwi_G0 Blackberry_G0       Kiwi_G0
## 6 Blackberry_Strawberry Blackberry_Strawberry_G0 Blackberry_G0 Strawberry_G0
##           pop_gen
## 1 Blackberry33_G0
## 2 Blackberry33_G0
## 3 Blackberry33_G0
## 4 Blackberry33_G0
## 5 Blackberry33_G0
## 6 Blackberry33_G0
levels(data_PREF$Population)
##  [1] "Blackberry31" "Blackberry32" "Blackberry33" "Blackberry34" "Blackberry35"
##  [6] "Blackberry36" "Blackberry37" "Blackberry38" "Blackberry39" "Blackberry40"
## [11] "Blackberry43" "Blackberry44" "Blackberry45" "Cherry103"    "Cherry104"   
## [16] "Cherry3"      "Cherry47"     "Cherry50"     "Cherry51"     "Cherry52"    
## [21] "Cherry6"      "Cherry7"      "Strawberry42" "Strawberry44" "Strawberry53"

2 OVIPOSITION STIMULATION

2.1 Exploration data

tapply(data_PERF$Nb_eggs, list(data_PERF$Original_environment, data_PERF$Test_environment, data_PERF$Generation), mean, na.rm = TRUE)
## , , G0
## 
##            Blackberry   Cherry Strawberry
## Blackberry    6.52459 10.22034   7.306122
## Cherry       23.35443 19.77500  12.283784
## Strawberry   21.30769 33.60000  16.833333
## 
## , , G2
## 
##            Blackberry   Cherry Strawberry
## Blackberry   103.9720 121.8491   98.71963
## Cherry       131.7174 133.8400  103.51064
## Strawberry   119.1739 119.4894  111.56522
ggplot2::ggplot(data = data_PERF[data_PERF$Generation=="G0",], 
                aes(x = Test_environment, y = Nb_eggs, color = Test_environment)) +
  facet_wrap( ~ Population) +
  geom_point() +
  geom_boxplot() +
  theme_LO_sober

ggplot2::ggplot(data = data_PERF, 
                aes(x = Nb_eggs, fill = Original_environment)) +
  facet_wrap( ~ Generation) +
  geom_histogram(position="identity", alpha=0.5) +
  theme_LO_sober

2.2 Analysis

#######################################################
## Analysis of genetic effects lm   ###
#######################################################
  
m1 <- aov(log(Nb_eggs+1) ~ pop_gen + hab_gen + SA + 
            Original_environment:Test_environment,
          contrasts = list(Original_environment = "contr.sum", 
                           Test_environment = "contr.sum"), 
          data = data_PERF)
  
  
## F test for SA
(Fratio <- (anova(m1)[3,2]/anova(m1)[4,2])/(1/anova(m1)[4, 1]))
## [1] 1.971666
(pvalue <- 1 - pf(Fratio, 1, anova(m1)[4, 1])) 
## [1] 0.2548912
#######################################################
## Analysis of non-genetic effects lm   ###
#######################################################
  
m2 <- aov(log(Nb_eggs+1) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0, 
          data = data_PERF)

## F test for SA
(Fratio_Gen <- (anova(m2)[3,2]/anova(m2)[5,2])/(1/anova(m2)[5, 1]))
## [1] 2.264412
(pvalue_Gen <- 1 - pf(Fratio_Gen, 1, anova(m2)[5, 1]) )
## [1] 0.2294357
## F test for SA
(Fratio_NonGen <- (anova(m2)[4,2]/anova(m2)[6,2])/(1/anova(m2)[6, 1]))
## [1] 1.369663
(pvalue_NonGen <- 1 - pf(Fratio_NonGen, 1, anova(m2)[6, 1]) )
## [1] 0.3263818
## Compute R2 for SA 
## Compute R2 = MS Interaction model without SA - MS Interaction model with SA / MS Interaction model without SA
(rsqgen <- 1-anova(m2)[5, 3]/((anova(m2)[3, 2]+anova(m2)[5, 2])/(anova(m2)[3, 1]+anova(m2)[5, 1])))
## [1] 0.240181
(rsqng <- 1-anova(m2)[6, 3]/((anova(m2)[4, 2]+anova(m2)[6, 2])/(anova(m2)[4, 1]+anova(m2)[6, 1])))
## [1] 0.08459751

##Plot

(PLOT_eggs_G0 <- plot_RTP_residuals(dataset = data_PERF, trait = "Nb_eggs", gen = "G0"))

(PLOT_eggs_G2 <- plot_RTP_residuals(dataset = data_PERF, trait = "Nb_eggs", gen = "G2"))

(PLOT_GEN_eggs_G0 <- plot_Genetic_Nongenetic_residuals(dataset = data_PERF, 
                                                   trait = "Nb_eggs", 
                                                   effect = "Non-genetic"))

(PLOT_GEN_eggs_G2 <- plot_Genetic_Nongenetic_residuals(dataset = data_PERF, 
                                                       trait = "Nb_eggs", 
                                                   effect = "Genetic"))

3 NUMBER OF ADULTS

3.1 Exploration data

tapply(data_PERF$Nb_adults, list(data_PERF$Original_environment, data_PERF$Test_environment, data_PERF$Generation), mean, na.rm = TRUE)
## , , G0
## 
##            Blackberry    Cherry Strawberry
## Blackberry   4.508197  4.457627   3.816327
## Cherry       6.873418  5.675000   2.378378
## Strawberry  12.230769 16.000000   5.416667
## 
## , , G2
## 
##            Blackberry   Cherry Strawberry
## Blackberry   27.31776 23.18868   13.15888
## Cherry       22.91304 23.32000   11.53191
## Strawberry   25.65217 22.36170   16.04348
ggplot2::ggplot(data = data_PERF[data_PERF$Generation=="G0",], 
                aes(x = Test_environment, y = Nb_adults, color = Test_environment)) +
  facet_wrap( ~ Population) +
  geom_point() +
  geom_boxplot() +
  theme_LO_sober

ggplot2::ggplot(data = data_PERF[data_PERF$Generation=="G2",], 
                aes(x = Test_environment, y = Nb_adults, color = Test_environment)) +
  facet_wrap( ~ Population) +
  geom_point() +
  geom_boxplot() +
  theme_LO_sober

ggplot2::ggplot(data = data_PERF, 
                aes(x = Nb_adults, fill = Original_environment)) +
  facet_wrap( ~ Generation) +
  geom_histogram(position="identity", alpha=0.5) +
  theme_LO_sober

3.2 Analysis

#######################################################
## Analysis of genetic effects lm   ###
#######################################################
  
m1 <- aov(log(Nb_adults+1) ~ pop_gen + hab_gen + SA + 
            Original_environment:Test_environment,
          contrasts = list(Original_environment = "contr.sum", 
                           Test_environment = "contr.sum"), 
          data = data_PERF)
  
  
## F test for SA
(Fratio <- (anova(m1)[3,2]/anova(m1)[4,2])/(1/anova(m1)[4, 1]))
## [1] 2.01345
(pvalue <- 1 - pf(Fratio, 1, anova(m1)[4, 1])) 
## [1] 0.2509626
#######################################################
## Analysis of non-genetic effects lm   ###
#######################################################
  
m2 <- aov(log(Nb_adults+1) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0, 
          data = data_PERF)

## F test for SA
(Fratio_Gen <- (anova(m2)[3,2]/anova(m2)[5,2])/(1/anova(m2)[5, 1]))
## [1] 2.447544
(pvalue_Gen <- 1 - pf(Fratio_Gen, 1, anova(m2)[5, 1]) )
## [1] 0.2156678
## F test for SA
(Fratio_NonGen <- (anova(m2)[4,2]/anova(m2)[6,2])/(1/anova(m2)[6, 1]))
## [1] 1.111569
(pvalue_NonGen <- 1 - pf(Fratio_NonGen, 1, anova(m2)[6, 1]) )
## [1] 0.3691494
## Compute R2 for SA 
## = MS Interaction model without SA - MS Interaction model with SA / MS Interaction model without SA
(r2_SA_genet <- 1-(anova(m2)[5, 3]/((anova(m2)[3, 2]+anova(m2)[5, 2])/(anova(m2)[3, 1]+anova(m2)[5, 1]))))
## [1] 0.2657241
(r2_SA_nongenet <- 1-(anova(m2)[6, 3]/((anova(m2)[4, 2]+anova(m2)[6, 2])/(anova(m2)[4, 1]+anova(m2)[6, 1]))))
## [1] 0.02713541

3.3 Analysis with eggs as covariable

#######################################################
## Should we consider the number of eggs?   ###
#######################################################

# Original
m1 <- aov(log(Nb_adults+1) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0, 
          data = data_PERF)


## Correction for number of eggs
m2 <- aov(log(Nb_adults+1) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 +
            log(Nb_eggs+1), 
          data = data_PERF)


## With egg score
m3 <- aov(log(Nb_adults+1) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 +
            EggScore, 
          data = data_PERF)

## Compare with 5 egg scores
m4 <- aov(log(Nb_adults+1) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 +
            EggScoreFive, 
          data = data_PERF)

## Compare with EggScoreSmall
m5 <- aov(log(Nb_adults+1) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 +
            EggScoreSmall, 
          data = data_PERF)



MuMIn::model.sel(m1, m2, m3, m4, m5)
## Model selection table 
##       (Int) hab_gen pop_gen SA IG0:SA Org_env:Tst_env IG0:Org_env:Tst_env
## m2 -0.06648       +       +  +      +               +                   +
## m5  0.41840       +       +  +      +               +                   +
## m3  0.42840       +       +  +      +               +                   +
## m4  0.44140       +       +  +      +               +                   +
## m1  0.62560       +       +  +      +               +                   +
##    log(Nb_egg+1) EgS ESF ESS             family df    logLik   AICc  delta
## m2        0.5125             gaussian(identity) 63  -977.876 2090.0   0.00
## m5                         + gaussian(identity) 67 -1040.074 2223.5 133.51
## m3                 +         gaussian(identity) 65 -1045.327 2229.5 139.45
## m4                     +     gaussian(identity) 66 -1044.861 2230.8 140.80
## m1                           gaussian(identity) 62 -1118.184 2368.4 278.35
##    weight
## m2      1
## m5      0
## m3      0
## m4      0
## m1      0
## Models ranked by AICc(x)
# Models are not all fitted to the same data: because 6 tubes without Nb_eggs are missing for m2, m3, m4 and m5


## Cl= The best model is when the number of eggs is considered as a continuous variable
### model m2 provides a better description of the data than model m1 






#######################################################
## Analysis of genetic effects lm   ###
#######################################################
  
m1 <- aov(log(Nb_adults+1) ~ pop_gen + hab_gen + SA + 
            Original_environment:Test_environment + log(Nb_eggs+1),
          contrasts = list(Original_environment = "contr.sum", 
                           Test_environment = "contr.sum"), 
          data = data_PERF)
  
  
## F test for SA
(Fratio <- (anova(m1)[3,2]/anova(m1)[5,2])/(1/anova(m1)[5, 1]))
## [1] 12.95037
(pvalue <- 1 - pf(Fratio, 1, anova(m1)[5, 1])) 
## [1] 0.03679697
#######################################################
## Analysis of non-genetic effects lm   ###
#######################################################
  
m2 <- aov(log(Nb_adults+1) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 + 
            log(Nb_eggs+1), 
          data = data_PERF)

## F test for SA
(Fratio_Gen <- (anova(m2)[3,2]/anova(m2)[6,2])/(1/anova(m2)[6, 1]))
## [1] 17.09331
(pvalue_Gen <- 1 - pf(Fratio_Gen, 1, anova(m2)[6, 1]) )
## [1] 0.02567856
## F test for SA
(Fratio_NonGen <- (anova(m2)[5,2]/anova(m2)[7,2])/(1/anova(m2)[7, 1]))
## [1] 0.6257893
(pvalue_NonGen <- 1 - pf(Fratio_NonGen, 1, anova(m2)[7, 1]) )
## [1] 0.4866763
## Compute R2 for SA 
## = MS Interaction model without SA - MS Interaction model with SA / MS Interaction model without SA
(r2_SA_genet <- 1-(anova(m2)[6, 3]/((anova(m2)[3, 2]+anova(m2)[6, 2])/(anova(m2)[3, 1]+anova(m2)[6, 1]))))
## [1] 0.8009287
(r2_SA_nongenet <- 1-(anova(m2)[7, 3]/((anova(m2)[5, 2]+anova(m2)[7, 2])/(anova(m2)[5, 1]+anova(m2)[7, 1]))))
## [1] -0.1032081

##Plot

(PLOT_adult_G0 <- plot_RTP_residuals(dataset = data_PERF, trait = "Nb_adults", gen = "G0"))

(PLOT_adult_G2 <- plot_RTP_residuals(dataset = data_PERF, trait = "Nb_adults", gen = "G2"))

(PLOT_GEN_adult_G0 <- plot_Genetic_Nongenetic_residuals(dataset = data_PERF, 
                                                    trait = "Nb_adults", 
                                                    effect = "Non-genetic"))

(PLOT_GEN_adult_G2 <- plot_Genetic_Nongenetic_residuals(dataset = data_PERF, 
                                                    trait = "Nb_adults", 
                                                    effect = "Genetic"))

4 OFFPSRING PERFORMANCE

4.1 Exploration data

tapply(data_PERF_Rate$Rate, list(data_PERF_Rate$Original_environment, 
                                 data_PERF_Rate$Test_environment, 
                                 data_PERF_Rate$Generation), mean, na.rm = TRUE)
## , , G0
## 
##            Blackberry    Cherry Strawberry
## Blackberry  0.6125868 0.4703813  0.5248347
## Cherry      0.2420863 0.3199467  0.2176492
## Strawberry  0.2888558 0.3595403  0.3373855
## 
## , , G2
## 
##            Blackberry    Cherry Strawberry
## Blackberry  0.2579090 0.1923140  0.1345134
## Cherry      0.1855892 0.1801664  0.1103531
## Strawberry  0.2206513 0.1950778  0.1619433
ggplot2::ggplot(data = data_PERF_Rate[data_PERF_Rate$Generation=="G0",], 
                aes(x = Test_environment, y = Rate, color = Test_environment)) +
  facet_wrap( ~ Population) +
  geom_point() +
  geom_boxplot() +
  theme_LO_sober

ggplot2::ggplot(data = data_PERF_Rate[data_PERF_Rate$Generation=="G2",], 
                aes(x = Test_environment, y = Rate, color = Test_environment)) +
  facet_wrap( ~ Population) +
  geom_point() +
  geom_boxplot() +
  theme_LO_sober

ggplot2::ggplot(data = data_PERF_Rate, 
                aes(x = Rate, fill = Original_environment)) +
  facet_wrap( ~ Generation) +
  geom_histogram(position="identity", alpha=0.5) +
  theme_LO_sober

lattice::xyplot(Rate~Nb_eggs|Original_environment*Test_environment, 
                data=data_PERF_Rate)

lattice::xyplot(Rate~EggScore|Original_environment*Test_environment, 
                data=data_PERF_Rate)

lattice::bwplot(Rate~EggScore|Original_environment*Test_environment, 
                data=data_PERF_Rate)

lattice::bwplot(Rate~EggScoreFive|Original_environment*Test_environment, 
       data=data_PERF_Rate)

lattice::bwplot(Rate~EggScoreSmall|Original_environment*Test_environment, 
       data=data_PERF_Rate)

AvEmergenceRate <- tapply(data_PERF_Rate$Rate, 
                          list(data_PERF_Rate$EggScoreFive,
                              data_PERF_Rate$Original_environment,
                              data_PERF_Rate$Test_environment),mean)
tapply(data_PERF_Rate$Rate, list(data_PERF_Rate$EggScoreFive,
                                 data_PERF_Rate$Original_environment,
                                 data_PERF_Rate$Test_environment), length)
## , , Blackberry
## 
##   Blackberry Cherry Strawberry
## 1         58     69         13
## 2         42     22         15
## 3         41     17         20
## 4         10     10          8
## 5          3      4          2
## 
## , , Cherry
## 
##   Blackberry Cherry Strawberry
## 1         45     70          8
## 2         31     18         19
## 3         50     21         21
## 4         23     11          9
## 5          2      7         NA
## 
## , , Strawberry
## 
##   Blackberry Cherry Strawberry
## 1         42     72         18
## 2         55     24         13
## 3         39     16         16
## 4          7      6          9
## 5         NA     NA          2
AvEmergenceRate[, , "Blackberry"][, "Cherry"]/AvEmergenceRate[, , "Blackberry"][, "Blackberry"]
##         1         2         3         4         5 
## 0.4133863 1.0258692 0.7334304 0.9735309 0.1709919
AvEmergenceRate[, , "Cherry"][, "Blackberry"]/AvEmergenceRate[, , "Cherry"][, "Cherry"]
##         1         2         3         4         5 
## 1.4425369 0.8562518 1.0861981 0.8514873 1.4918117
AvEmergenceRate[, , "Blackberry"][, "Cherry"]/AvEmergenceRate[, , "Cherry"][, "Cherry"]
##         1         2         3         4         5 
## 0.7219048 1.0522488 1.0774054 0.9037047 0.4651315

4.2 Analysis

#######################################################
## Analysis of genetic effects lm   ###
#######################################################
  
m1 <- aov(asin(sqrt(Rate)) ~ pop_gen + hab_gen + SA + 
            Original_environment:Test_environment,
          contrasts = list(Original_environment = "contr.sum", 
                           Test_environment = "contr.sum"), 
          data = data_PERF_Rate)
  
  
## F test for SA
(Fratio <- (anova(m1)[3,2]/anova(m1)[4,2])/(1/anova(m1)[4, 1]))
## [1] 30.89092
(pvalue <- 1 - pf(Fratio, 1, anova(m1)[4, 1])) 
## [1] 0.0114893
#######################################################
## Analysis of non-genetic effects lm   ###
#######################################################
  
m2 <- aov(asin(sqrt(Rate)) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0, 
          data = data_PERF_Rate)

## F test for SA
(Fratio_Gen <- (anova(m2)[3,2]/anova(m2)[5,2])/(1/anova(m2)[5, 1]))
## [1] 59.30375
(pvalue_Gen <- 1 - pf(Fratio_Gen, 1, anova(m2)[5, 1]) )
## [1] 0.004550853
## F test for SA
(Fratio_NonGen <- (anova(m2)[4,2]/anova(m2)[6,2])/(1/anova(m2)[6, 1]))
## [1] 3.305405
(pvalue_NonGen <- 1 - pf(Fratio_NonGen, 1, anova(m2)[6, 1]) )
## [1] 0.1666402
## Compute R2 for SA 
## = MS Interaction model without SA - MS Interaction model with SA / MS Interaction model without SA
(r2_SA_genet <- 1-(anova(m2)[5, 3]/((anova(m2)[3, 2]+anova(m2)[5, 2])/(anova(m2)[3, 1]+anova(m2)[5, 1]))))
## [1] 0.9357984
(r2_SA_nongenet <- 1-(anova(m2)[6, 3]/((anova(m2)[4, 2]+anova(m2)[6, 2])/(anova(m2)[4, 1]+anova(m2)[6, 1]))))
## [1] 0.3656236

4.3 Analysis with eggs as covariable

#######################################################
## Should we consider the number of eggs?   ###
#######################################################

# Original
m1 <- aov(asin(sqrt(Rate)) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0, 
          data = data_PERF_Rate)



## Correction for number of eggs
m2 <- aov(asin(sqrt(Rate)) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 +
            log(Nb_eggs), 
          data = data_PERF_Rate)




## With egg score
m3 <- aov(asin(sqrt(Rate)) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 +
            EggScore, 
          data = data_PERF_Rate)

## Compare with 5 egg scores
m4 <- aov(asin(sqrt(Rate)) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 +
            EggScoreFive, 
          data = data_PERF_Rate)

## Compare with EggScoreSmall
m5 <- aov(asin(sqrt(Rate)) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 +
            EggScoreSmall, 
          data = data_PERF_Rate)

## Correction for number of eggs
m6 <- aov(asin(sqrt(Rate)) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 +
            Nb_eggs, 
          data = data_PERF_Rate)


data_PERF_Rate$Square_NbEggs <- data_PERF_Rate$Nb_eggs^2
## Correction for number of eggs
m7 <- aov(asin(sqrt(Rate)) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 +
            Nb_eggs + Square_NbEggs, 
          data = data_PERF_Rate)


## Correction for number of eggs
m8 <- aov(asin(sqrt(Rate)) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 +
            Square_NbEggs, 
          data = data_PERF_Rate)

MuMIn::model.sel(m1, m2, m3, m4, m5, m6, m7, m8)
## Model selection table 
##     (Int) hab_gen pop_gen SA IG0:SA Org_env:Tst_env IG0:Org_env:Tst_env
## m2 0.6870       +       +  +      +               +                   +
## m7 0.6447       +       +  +      +               +                   +
## m6 0.6187       +       +  +      +               +                   +
## m8 0.6060       +       +  +      +               +                   +
## m1 0.6144       +       +  +      +               +                   +
## m3 0.6155       +       +  +      +               +                   +
## m4 0.6145       +       +  +      +               +                   +
## m5 0.6159       +       +  +      +               +                   +
##    log(Nb_egg) EgS ESF ESS    Nb_egg    Sqr_NbE             family df   logLik
## m2    -0.06172                                  gaussian(identity) 63 -135.321
## m7                         -0.002675  7.317e-06 gaussian(identity) 64 -141.008
## m6                         -0.001058            gaussian(identity) 63 -143.208
## m8                                   -3.155e-06 gaussian(identity) 63 -146.200
## m1                                              gaussian(identity) 62 -149.207
## m3               +                              gaussian(identity) 65 -146.871
## m4                   +                          gaussian(identity) 66 -146.858
## m5                       +                      gaussian(identity) 67 -146.342
##     AICc delta weight
## m2 405.4  0.00  0.999
## m7 419.0 13.66  0.001
## m6 421.1 15.77  0.000
## m8 427.1 21.76  0.000
## m1 430.9 25.49  0.000
## m3 433.0 27.68  0.000
## m4 435.3 29.95  0.000
## m5 436.6 31.22  0.000
## Models ranked by AICc(x)
## Cl= The best model is when the number of eggs is considered as a continuous variable
### model m2 provides a better description of the data than model m1 







#######################################################
## Analysis of genetic effects lm   ###
#######################################################
  
m1 <- aov(asin(sqrt(Rate)) ~ pop_gen + hab_gen + SA + 
            Original_environment:Test_environment +
            log(Nb_eggs),
          contrasts = list(Original_environment = "contr.sum", 
                           Test_environment = "contr.sum"), 
          data = data_PERF_Rate)
  
  
## F test for SA
(Fratio <- (anova(m1)[3,2]/anova(m1)[5,2])/(1/anova(m1)[5, 1]))
## [1] 27.56268
(pvalue <- 1 - pf(Fratio, 1, anova(m1)[5, 1])) 
## [1] 0.01345815
## Compute R2 = MS Interaction model without SA - MS Interaction model with SA / MS Interaction model without SA
(rsqgen <- 1-anova(m1)[5, 3]/((anova(m1)[3, 2]+anova(m1)[4, 2])/(anova(m1)[3, 1]+anova(m1)[4, 1])))
## [1] 0.9762833
#######################################################
## Analysis of non-genetic effects lm   ###
#######################################################
  
m2 <- aov(asin(sqrt(Rate)) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 +
            log(Nb_eggs), 
          data = data_PERF_Rate)


## F test for SA
(Fratio_Gen <- (anova(m2)[3,2]/anova(m2)[6,2])/(1/anova(m2)[6, 1]))
## [1] 37.67609
(pvalue_Gen <- 1 - pf(Fratio_Gen, 1, anova(m2)[6, 1]) )
## [1] 0.008696739
## F test for SA
(Fratio_NonGen <- (anova(m2)[5,2]/anova(m2)[7,2])/(1/anova(m2)[7, 1]))
## [1] 2.164733
(pvalue_NonGen <- 1 - pf(Fratio_NonGen, 1, anova(m2)[7, 1]) )
## [1] 0.2375861
# Compute R2 = MS Interaction model without SA - MS Interaction model with SA / MS Interaction model without SA
rsqgen <- 1-anova(m2)[6, 3]/((anova(m2)[3, 2]+anova(m2)[6, 2])/(anova(m2)[3, 1]+anova(m2)[6, 1]))
rsqng <- 1-anova(m2)[7, 3]/((anova(m2)[5, 2]+anova(m2)[7, 2])/(anova(m2)[5, 1]+anova(m2)[7, 1]))
  

## Compute R2 for SA 
## = MS Interaction model without SA - MS Interaction model with SA / MS Interaction model without SA
(r2_SA_genet <- 1-(anova(m2)[6, 3]/((anova(m2)[3, 2]+anova(m2)[6, 2])/(anova(m2)[3, 1]+anova(m2)[6, 1]))))
## [1] 0.9016621
(r2_SA_nongenet <- 1-(anova(m2)[7, 3]/((anova(m2)[5, 2]+anova(m2)[7, 2])/(anova(m2)[5, 1]+anova(m2)[7, 1]))))
## [1] 0.2255167

4.4 Heterogeneity across populations

data_Rate_G2 <- data_PERF_Rate[data_PERF_Rate$Generation=="G2",]
data_Rate_G2 <- data_Rate_G2[complete.cases(data_Rate_G2$Rate), ]

  
m0 <- lm(Rate~Original_environment*Test_environment, data=data_Rate_G2,)

summary(m0)
## 
## Call:
## lm(formula = Rate ~ Original_environment * Test_environment, 
##     data = data_Rate_G2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.25791 -0.09879 -0.02409  0.07267  0.59923 
## 
## Coefficients:
##                                                           Estimate Std. Error
## (Intercept)                                                0.25791    0.01305
## Original_environmentCherry                                -0.07232    0.02372
## Original_environmentStrawberry                            -0.03726    0.02372
## Test_environmentCherry                                    -0.06560    0.01845
## Test_environmentStrawberry                                -0.12340    0.01841
## Original_environmentCherry:Test_environmentCherry          0.06017    0.03307
## Original_environmentStrawberry:Test_environmentCherry      0.04002    0.03342
## Original_environmentCherry:Test_environmentStrawberry      0.04816    0.03339
## Original_environmentStrawberry:Test_environmentStrawberry  0.06469    0.03352
##                                                           t value Pr(>|t|)    
## (Intercept)                                                19.767  < 2e-16 ***
## Original_environmentCherry                                 -3.049 0.002397 ** 
## Original_environmentStrawberry                             -1.571 0.116740    
## Test_environmentCherry                                     -3.555 0.000408 ***
## Test_environmentStrawberry                                 -6.703 4.76e-11 ***
## Original_environmentCherry:Test_environmentCherry           1.820 0.069337 .  
## Original_environmentStrawberry:Test_environmentCherry       1.198 0.231534    
## Original_environmentCherry:Test_environmentStrawberry       1.442 0.149775    
## Original_environmentStrawberry:Test_environmentStrawberry   1.930 0.054091 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1343 on 592 degrees of freedom
## Multiple R-squared:  0.1004, Adjusted R-squared:  0.08825 
## F-statistic: 8.259 on 8 and 592 DF,  p-value: 1.223e-10
tapply(data_Rate_G2$Rate, 
       list(data_Rate_G2$Original_environment, 
            data_Rate_G2$Test_environment), mean)
##            Blackberry    Cherry Strawberry
## Blackberry  0.2579090 0.1923140  0.1345134
## Cherry      0.1855892 0.1801664  0.1103531
## Strawberry  0.2206513 0.1950778  0.1619433
tapply(data_Rate_G2$Rate, 
       list(data_Rate_G2$Original_environment, 
            data_Rate_G2$Test_environment), var)
##            Blackberry      Cherry Strawberry
## Blackberry 0.03005193 0.013749525 0.01436498
## Cherry     0.02922031 0.019467766 0.01399160
## Strawberry 0.01555929 0.009464032 0.01139614
range(data_Rate_G2$Nb_adults, na.rm = TRUE)
## [1]   0 108
## Check number of eggs and adults
tapply(data_Rate_G2$Nb_eggs, 
       list(data_Rate_G2$Original_environment, 
            data_Rate_G2$Test_environment), mean)
##            Blackberry   Cherry Strawberry
## Blackberry   104.2170 121.8491   98.71963
## Cherry       131.7174 133.8400  103.51064
## Strawberry   119.1739 119.4894  111.56522
tapply(data_Rate_G2$Nb_adults, 
       list(data_Rate_G2$Original_environment, 
            data_Rate_G2$Test_environment), mean)
##            Blackberry   Cherry Strawberry
## Blackberry   26.82075 23.18868   13.15888
## Cherry       22.91304 23.32000   11.53191
## Strawberry   25.65217 22.36170   16.04348
## Check for the presence of negative correlations
m1 <- lm(asin(sqrt(Rate)) ~ Population + Test_environment,
         data=data_Rate_G2)
data_Rate_G2$resid <- residuals(m1)

meanbypopbytestenv <- as.data.frame(tapply(data_Rate_G2$resid, 
                                           list(data_Rate_G2$Population,
                                                data_Rate_G2$Test_environment), mean))


## Cherry ~ Blackberry
plot(meanbypopbytestenv$Blackberry, 
     meanbypopbytestenv$Cherry, 
     col=as.numeric(data_Rate_G2$Original_environment[match(rownames(meanbypopbytestenv),
                                                    data_Rate_G2$Population)]), 
     xlab="Blackberry", ylab="Cherry", pch=16)
legend("topright", levels(data_Rate_G2$Original_environment), 
       col=as.numeric(as.factor(levels(data_Rate_G2$Original_environment))), pch=16)

## Strawberry ~ Cherry
plot(meanbypopbytestenv$Cherry,  meanbypopbytestenv$Strawberry, col=as.numeric(data_Rate_G2$Original_environment[match(rownames(meanbypopbytestenv), data_Rate_G2$Population)]), xlab="Cherry", ylab="Strawberry", pch=16)
legend("bottomright", levels(data_Rate_G2$Original_environment), col=as.numeric(as.factor(levels(data_Rate_G2$Original_environment))), pch=16)

## Strawberry ~ Blackberry
plot(meanbypopbytestenv$Blackberry, meanbypopbytestenv$Strawberry, col=as.numeric(data_Rate_G2$Original_environment[match(rownames(meanbypopbytestenv), data_Rate_G2$Population)]), xlab="Blackberry", ylab="Strawberry", pch=16)
legend("topright", levels(data_Rate_G2$Original_environment), col=as.numeric(as.factor(levels(data_Rate_G2$Original_environment))), pch=16)

##Plot

(PLOT_rate_G0 <- plot_RTP_residuals(dataset = data_PERF_Rate, trait = "Rate", gen = "G0"))

(PLOT_rate_G2 <- plot_RTP_residuals(dataset = data_PERF_Rate, trait = "Rate", gen = "G2"))

(PLOT_GEN_rate_G0 <- plot_Genetic_Nongenetic_residuals(dataset = data_PERF_Rate, 
                                                    trait = "Rate", 
                                                    effect = "Non-genetic"))

(PLOT_GEN_rate_G2 <- plot_Genetic_Nongenetic_residuals(dataset = data_PERF_Rate, 
                                                    trait = "Rate", 
                                                    effect = "Genetic"))

5 OVIOSITION PREFERENCE

5.1 Twelve fruits

5.1.1 Exploration data

tapply(data_PREF$Nb_eggs, list(data_PREF$Original_environment, 
                                 data_PREF$Test_environment, 
                                 data_PREF$Generation), mean, na.rm = TRUE)
## , , G0
## 
##              Apricot Blackberry Blackcurrant    Cherry Cranberry       Fig
## Blackberry 0.2758621  0.7931034    0.1379310 0.4137931 0.2068966 0.2586207
## Cherry     0.3300000  0.5000000    1.2600000 2.4800000 0.3800000 0.7400000
## Strawberry 0.5714286  1.1904762    0.6666667 1.6666667 0.5714286 0.4285714
##                Grape      Kiwi Raspberry  Rosehips Strawberry    Tomato
## Blackberry 0.2068966 0.0862069 0.4655172 0.3793103  0.5862069 0.2241379
## Cherry     0.9700000 1.0500000 1.5800000 1.2400000  0.6700000 0.6969697
## Strawberry 0.4761905 0.1428571 0.6666667 0.3809524  0.8571429 0.0952381
## 
## , , G2
## 
##              Apricot Blackberry Blackcurrant   Cherry Cranberry      Fig
## Blackberry  6.826531   24.62245     18.61224 23.41837  8.428571 12.73469
## Cherry     10.423077   16.09615     19.09615 30.38462  7.480769 13.46154
## Strawberry 15.450000   21.92500     22.12500 24.92500  7.200000 13.02500
##               Grape     Kiwi Raspberry Rosehips Strawberry   Tomato
## Blackberry 24.55102 16.57143  16.59184 11.23469   6.938776 12.53061
## Cherry     11.67308 11.90385  12.78846 11.46154  10.673077 10.59615
## Strawberry 12.55000 14.87500  19.45000 16.20000  12.775000 14.90000
ggplot2::ggplot(data = data_PREF, 
                aes(x = Test_environment, y = Nb_eggs, color = Test_environment)) +
  facet_wrap( ~ Population) +
  geom_point() +
  geom_boxplot() +
  theme_LO_sober

ggplot2::ggplot(data = data_PREF, 
                aes(x = Nb_eggs, fill = Original_environment)) +
  facet_wrap( ~ Generation) +
  geom_histogram(position="identity", alpha=0.5) +
  theme_LO_sober

5.1.2 Analysis

#######################################################
## Analysis of genetic effects lm   ###
#######################################################
  
m1 <- aov(log(Nb_eggs+1) ~ pop_gen + hab_gen + SA + 
            Original_environment:Test_environment,
          contrasts = list(Original_environment = "contr.sum", 
                           Test_environment = "contr.sum"), 
          data = data_PREF)
  
  
## F test for SA
(Fratio <- (anova(m1)[3,2]/anova(m1)[4,2])/(1/anova(m1)[4, 1]))
## [1] 10.1742
(pvalue <- 1 - pf(Fratio, 1, anova(m1)[4, 1])) 
## [1] 0.004407544
#######################################################
## Analysis of non-genetic effects lm   ###
#######################################################
  
m2 <- aov(log(Nb_eggs+1) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 +
            BoxID, 
          data = data_PREF)

## F test for SA
(Fratio_Gen <- (anova(m2)[3,2]/anova(m2)[6,2])/(1/anova(m2)[6, 1]))
## [1] 9.833047
(pvalue_Gen <- 1 - pf(Fratio_Gen, 1, anova(m2)[6, 1]) )
## [1] 0.004993505
## F test for SA
(Fratio_NonGen <- (anova(m2)[5,2]/anova(m2)[7,2])/(1/anova(m2)[7, 1]))
## [1] 2.500166
(pvalue_NonGen <- 1 - pf(Fratio_NonGen, 1, anova(m2)[7, 1]) )
## [1] 0.1287796
## Compute R2 for SA 
## = MS Interaction model without SA - MS Interaction model with SA / MS Interaction model without SA
(r2_SA_genet <- 1-(anova(m2)[6, 3]/((anova(m2)[3, 2]+anova(m2)[6, 2])/(anova(m2)[3, 1]+anova(m2)[6, 1]))))
## [1] 0.2864799
(r2_SA_nongenet <- 1-(anova(m2)[7, 3]/((anova(m2)[5, 2]+anova(m2)[7, 2])/(anova(m2)[5, 1]+anova(m2)[7, 1]))))
## [1] 0.0638364
#Local adaptation pattern: 
lm_val = lm(log(Nb_eggs+1) ~ Test_environment + Population + SA + 
                      Test_environment:Original_environment + BoxID, 
                    data = data_PREF[data_PREF$Generation=="G0",])

(Fratio = anova(lm_val)[3,3]/anova(lm_val)[5,3])
## [1] 4.476081
(pvalue = 1 - pf(Fratio,anova(lm_val)[3,1],anova(lm_val)[5,1]))
## [1] 0.0464917
(df1 = anova(lm_val)[3,1])
## [1] 1
(df2 = anova(lm_val)[5,1])
## [1] 21
lm_val = lm(log(Nb_eggs+1) ~ Test_environment + Population + SA + 
                      Test_environment:Original_environment + BoxID, 
                    data = data_PREF[data_PREF$Generation=="G2",])

(Fratio = anova(lm_val)[3,3]/anova(lm_val)[5,3])
## [1] 6.319972
(pvalue = 1 - pf(Fratio,anova(lm_val)[3,1],anova(lm_val)[5,1]))
## [1] 0.02016054
(df1 = anova(lm_val)[3,1])
## [1] 1
(df2 = anova(lm_val)[5,1])
## [1] 21

5.2 Three fruits

5.2.1 Exploration data

tapply(data_PREF_three$Nb_eggs, list(data_PREF_three$Original_environment, 
                                 data_PREF_three$Test_environment, 
                                 data_PREF_three$Generation), mean, na.rm = TRUE)
## , , G0
## 
##            Blackberry    Cherry Strawberry
## Blackberry  0.7931034 0.4137931  0.5862069
## Cherry      0.5000000 2.4800000  0.6700000
## Strawberry  1.1904762 1.6666667  0.8571429
## 
## , , G2
## 
##            Blackberry   Cherry Strawberry
## Blackberry   24.62245 23.41837   6.938776
## Cherry       16.09615 30.38462  10.673077
## Strawberry   21.92500 24.92500  12.775000
ggplot2::ggplot(data = data_PREF_three, 
                aes(x = Test_environment, y = Nb_eggs, color = Test_environment)) +
  facet_wrap( ~ Population) +
  geom_point() +
  geom_boxplot() +
  theme_LO_sober

ggplot2::ggplot(data = data_PREF_three, 
                aes(x = Nb_eggs, fill = Original_environment)) +
  facet_wrap( ~ Generation) +
  geom_histogram(position="identity", alpha=0.5) +
  theme_LO_sober

5.2.2 Analysis

#######################################################
## Analysis of genetic effects lm   ###
#######################################################
  
m1 <- aov(log(Nb_eggs+1) ~ pop_gen + hab_gen + SA + 
            Original_environment:Test_environment +
            BoxID,
          contrasts = list(Original_environment = "contr.sum", 
                           Test_environment = "contr.sum"), 
          data = data_PREF_three)
  
  
## F test for SA
(Fratio <- (anova(m1)[3,2]/anova(m1)[5,2])/(1/anova(m1)[5, 1]))
## [1] 18.38385
(pvalue <- 1 - pf(Fratio, 1, anova(m1)[5, 1])) 
## [1] 0.02331823
#######################################################
## Analysis of non-genetic effects lm   ###
#######################################################
  
m2 <- aov(log(Nb_eggs+1) ~ pop_gen + hab_gen + SA:IndicG0 + SA +
            Original_environment:Test_environment + 
            Original_environment:Test_environment:IndicG0 +
            BoxID, 
          data = data_PREF_three)

## F test for SA
(Fratio_Gen <- (anova(m2)[3,2]/anova(m2)[6,2])/(1/anova(m2)[6, 1]))
## [1] 12.77418
(pvalue_Gen <- 1 - pf(Fratio_Gen, 1, anova(m2)[6, 1]) )
## [1] 0.0374429
## F test for SA
(Fratio_NonGen <- (anova(m2)[5,2]/anova(m2)[7,2])/(1/anova(m2)[7, 1]))
## [1] 3.536493
(pvalue_NonGen <- 1 - pf(Fratio_NonGen, 1, anova(m2)[7, 1]) )
## [1] 0.1566089
## Compute R2 for SA 
# ## = MS Interaction model without SA - MS Interaction model with SA / MS Interaction model without SA
# (r2_SA_genet <- 1-(anova(m2)[6, 3]/((anova(m2)[3, 2]+anova(m2)[6, 2])/(anova(m2)[3, 1]+anova(m2)[6, 1]))))
# 
# (r2_SA_nongenet <- 1-(anova(m2)[7, 3]/((anova(m2)[5, 2]+anova(m2)[7, 2])/(anova(m2)[5, 1]+anova(m2)[7, 1]))))
  ## Compute R2 = MS Interaction model without SA - MS Interaction model with SA / MS Interaction model without SA
rsqgen <- 1-anova(m2)[6, 3]/((anova(m2)[3, 2]+anova(m2)[6, 2])/(anova(m2)[3, 1]+anova(m2)[6, 1]))
rsqng <- 1-anova(m2)[7, 3]/((anova(m2)[5, 2]+anova(m2)[7, 2])/(anova(m2)[5, 1]+anova(m2)[7, 1]))
  rsqgen
## [1] 0.746421
  rsqng
## [1] 0.3880511

5.2.3 Plot

(PLOT_pref_G0 <- plot_RTP_residuals(dataset = data_PREF_three, 
                                    trait = "Nb_eggs", gen = "G0"))

(PLOT_pref_G2 <- plot_RTP_residuals(dataset = data_PREF_three, 
                                    trait = "Nb_eggs", gen = "G2"))

(PLOT_GEN_pref_G0 <- plot_Genetic_Nongenetic_residuals(dataset = data_PREF_three, 
                                    trait = "Nb_eggs", effect = "Non-genetic"))

(PLOT_GEN_pref_G2 <- plot_Genetic_Nongenetic_residuals(dataset = data_PREF_three, 
                                    trait = "Nb_eggs", effect = "Genetic"))

6 PLOT LOCAL ADAPTATION

legend <- lemon::g_legend(PLOT_pref_G0)
# 
# ############## FIRST GENERATION 
# LOCAL_ADAPTATION_PLOT_FRST <- cowplot::ggdraw() +
#   cowplot::draw_plot(PLOT_eggs_G0+theme(legend.position = "none", 
#                                         plot.title = element_blank()),
#             x =0, y = 0.5, width = 0.36, height = 0.46) +
#   cowplot::draw_plot(PLOT_adult_G0+theme(legend.position = "none", 
#                                         plot.title = element_blank()), 
#             x = 0.4, y = 0.5, width =0.36, height = 0.46) +
#   cowplot::draw_plot(PLOT_rate_G0+theme(legend.position = "none", 
#                                         plot.title = element_blank()),
#             x =0, y = 0, width = 0.36, height = 0.46) +
#   cowplot::draw_plot(PLOT_pref_G0+theme(legend.position = "none", 
#                                         plot.title = element_blank()), 
#             x = 0.4, y = 0, width = 0.36, height = 0.46) +
#   cowplot::draw_plot(legend, x = 0.85, y = 0.5, width = 0.0001, height = 0.0001) +
#   cowplot::draw_plot_label(c("A", "B", "C", "D"), 
#                            c(0.01, 0.4,0.01, 0.4), 
#                   c(1, 1,0.5,0.5), size = 19)
# # + cowplot::draw_plot_label("First generation", x = 0.2, y = 1 , size = 14) 
# LOCAL_ADAPTATION_PLOT_FRST
# 
# 
# # cowplot::save_plot(file =here::here("figures","Reciprocal_experiment_FirstGeneration.pdf"),
# #                   LOCAL_ADAPTATION_PLOT_FRST, base_height = 20/cm(1),
# #                   base_width = 30/cm(1), dpi = 1200)
# 
# 
# ############## THIRD GENERATION 
# LOCAL_ADAPTATION_PLOT_THIRD <- cowplot::ggdraw() +
#     cowplot::draw_plot(PLOT_eggs_G0+theme(legend.position = "none", 
#                                         plot.title = element_blank()),
#             x =0, y = 0.5, width = 0.36, height = 0.46) +
#   cowplot::draw_plot(PLOT_rate_G0+theme(legend.position = "none", 
#                                         plot.title = element_blank()), 
#             x = 0.4, y = 0.5, width =0.36, height = 0.46) +
#   cowplot::draw_plot(PLOT_rate_G0+theme(legend.position = "none", 
#                                         plot.title = element_blank()),
#             x =0, y = 0, width = 0.36, height = 0.46) +
#   cowplot::draw_plot(PLOT_pref_G0+theme(legend.position = "none", 
#                                         plot.title = element_blank()),
#             x =0, y = 0.5, width = 0.36, height = 0.46) +
#   cowplot::draw_plot(PLOT_adult_G2+theme(legend.position = "none", 
#                                         plot.title = element_blank()), 
#             x = 0.4, y = 0.5, width =0.36, height = 0.46) +
#   cowplot::draw_plot(PLOT_rate_G2+theme(legend.position = "none", 
#                                         plot.title = element_blank()),
#             x =0, y = 0, width = 0.36, height = 0.46) +
#   cowplot::draw_plot(PLOT_pref_G2+theme(legend.position = "none", 
#                                         plot.title = element_blank()), 
#             x = 0.4, y = 0, width = 0.36, height = 0.46) +
#   cowplot::draw_plot(legend, x = 0.85, y = 0.5, width = 0.0001, height = 0.0001) +
#   cowplot::draw_plot_label(c("A", "B", "C", "D"), 
#                            c(0.01, 0.4,0.01, 0.4), 
#                   c(1, 1,0.5,0.5), size = 19)
# LOCAL_ADAPTATION_PLOT_THIRD
# # 
# 
# cowplot::save_plot(file = here::here("figures","Reciprocal_experiment_ThirdGeneration.pdf"),
#                   LOCAL_ADAPTATION_PLOT_THIRD, base_height = 20/cm(1),
#                   base_width = 30/cm(1), dpi = 1200)
# 



## ALL GENERATIONS 
LOCAL_ADAPTATION_PLOT <- cowplot::ggdraw() + 
  cowplot::draw_plot(PLOT_pref_G0 + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.01, y = 0.5, width = 0.25, height = 0.45) + 
  cowplot::draw_plot(PLOT_eggs_G0 + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.31, y = 0.5, width = 0.25, height = 0.45) + 
  cowplot::draw_plot(PLOT_rate_G0 + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.61, y = 0.5, width = 0.25, height = 0.45) + 
  cowplot::draw_plot(legend, x = 0.93, y = 0.5, width = 0.0001, height = 0.0001) + 
  cowplot::draw_plot(PLOT_pref_G2 + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.01, y = 0, width = 0.25, height = 0.45) + 
  cowplot::draw_plot(PLOT_eggs_G2 + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.31, y = 0, width = 0.25, height = 0.45) + 
  cowplot::draw_plot(PLOT_rate_G2 + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.61, y = 0, width = 0.25, height = 0.45) + 
  cowplot::draw_plot_label(c("Generation: G0","Generation: G0","Generation: G1", "A", "B", "C", " ",
                    "Generation: G2","Generation: G2","Generation: G3", "D", "E", "F", " "),  
                  x = c(0.1,0.4,0.7, 0.01, 0.30, 0.61, 0.92, 0.10,0.4,0.7, 0.01, 0.3, 0.61, 0.92), 
                  y = c(1,1,1, 0.98, 0.98, 0.98, 0.98, 0.5, 0.5, 0.5, 0.48, 0.48, 0.48, 0.48), 
                  hjust = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0), 
                  vjust = c(1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5),
                  size = 16) 
 LOCAL_ADAPTATION_PLOT

# 
#  cowplot::save_plot(file =here::here("figures", "LOCAL_ADAPTATION_G0_G2.pdf"),
#                    LOCAL_ADAPTATION_PLOT,
#                    base_height = 18/cm(1), base_width = 34/cm(1), dpi = 610)
# 
# 

7 PLOT GENETIC vs PLASTIC

legend <- lemon::g_legend(PLOT_pref_G0)

Genetic_NonGenetic_PLOT <- cowplot::ggdraw() + 
  cowplot::draw_plot(PLOT_GEN_pref_G2 + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.01, y = 0.5, width = 0.25, height = 0.45) + 
  cowplot::draw_plot(PLOT_GEN_eggs_G2 + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.31, y = 0.5, width = 0.25, height = 0.45) + 
  cowplot::draw_plot(PLOT_GEN_rate_G2 + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.61, y = 0.5, width = 0.25, height = 0.45) + 
  cowplot::draw_plot(legend, x = 0.93, y = 0.5, width = 0.0001, height = 0.0001) + 
  cowplot::draw_plot(PLOT_GEN_pref_G0 + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.01, y = 0, width = 0.25, height = 0.45) + 
  cowplot::draw_plot(PLOT_GEN_eggs_G0 + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.31, y = 0, width = 0.25, height = 0.45) + 
  cowplot::draw_plot(PLOT_GEN_rate_G0 + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.61, y = 0, width = 0.25, height = 0.45) + 
  cowplot::draw_plot_label(c("Genetic effects", "A", "B", "C", " ",
                    "Plastic effects", "D", "E", "F", " "),  
                  x = c(0.4, 0.01, 0.30, 0.61, 0.92, 0.4, 0.01, 0.3, 0.61, 0.92), 
                  y = c(1, 0.98, 0.98, 0.98, 0.98, 0.5, 0.48, 0.48, 0.48, 0.48), 
                  hjust = c(0,0,0,0,0,0,0,0,0,0), 
                  vjust = c(1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5),
                  size = 16) 
 Genetic_NonGenetic_PLOT

# 
#
# cowplot::save_plot(file =here::here("figures", "GENETIC_NONGENETIC_First_Third.pdf"),
#                   Genetic_NonGenetic_PLOT,
#                   base_height = 18/cm(1), base_width = 34/cm(1), dpi = 610)

#

8 Plot old heterogeneity

################# STIMULATION 
###PLOT 
# (PAIR_BLACK_CHERRY_G2_EGGS <- plot_PairwisePOP_residuals(dataset = data_PERF, 
#                                                            trait = "Nb_eggs", gen = "G2", 
#                                                            fruit1 = "Cherry", fruit2 = "Blackberry"))
# 
# (PAIR_CHERRY_STRAW_G2_EGGS  <- plot_PairwisePOP_residuals(dataset = data_PERF, 
#                                                             trait = "Nb_eggs", gen = "G2",
#                                                             fruit1 = "Strawberry", fruit2 = "Cherry"))
# 
# (PAIR_STRW_BLACK_G2_EGGS  <- plot_PairwisePOP_residuals(dataset = data_PERF, 
#                                                           trait = "Nb_eggs", gen = "G2",
#                                                           fruit1 = "Blackberry", fruit2 = "Strawberry"))
# 
# 
# 
# 
# ###PLOT 
# (PAIR_BLACK_CHERRY_EGGS <- plot_PairwisePOP_residuals(dataset = data_PERF, 
#                                                            trait = "Nb_eggs", gen = "Both", 
#                                                            fruit1 = "Cherry", fruit2 = "Blackberry"))
# 
# (PAIR_CHERRY_STRAW_EGGS  <- plot_PairwisePOP_residuals(dataset = data_PERF, 
#                                                             trait = "Nb_eggs", gen = "Both",
#                                                             fruit1 = "Strawberry", fruit2 = "Cherry"))
# 
# (PAIR_STRW_BLACK_EGGS  <- plot_PairwisePOP_residuals(dataset = data_PERF, 
#                                                           trait = "Nb_eggs", gen = "Both",
#                                                           fruit1 = "Blackberry", 
#                                                      fruit2 = "Strawberry"))
# 
# 
# ################# PERFORMANCE 
# ###PLOT 
# (PAIR_BLACK_CHERRY_G2_RATE <- plot_PairwisePOP_residuals(dataset = data_PERF_Rate, trait = "Rate", gen = "G2", 
#                                        fruit1 = "Cherry", fruit2 = "Blackberry"))
# 
# (PAIR_CHERRY_STRAW_G2_RATE <- plot_PairwisePOP_residuals(dataset = data_PERF_Rate, trait = "Rate", gen = "G2", 
#                                        fruit1 = "Strawberry", fruit2 = "Cherry"))
# 
# (PAIR_STRW_BLACK_G2_RATE <- plot_PairwisePOP_residuals(dataset = data_PERF_Rate, trait = "Rate", gen = "G2", 
#                                        fruit1 = "Blackberry", fruit2 = "Strawberry"))
# 
# 
# (PAIR_BLACK_CHERRY_RATE <- plot_PairwisePOP_residuals(dataset = data_PERF_Rate, 
#                                                       trait = "Rate", 
#                                                       gen = "Both", 
#                                                       fruit1 = "Cherry", 
#                                                       fruit2 = "Blackberry"))
# 
# (PAIR_CHERRY_STRAW_RATE <- plot_PairwisePOP_residuals(dataset = data_PERF_Rate, 
#                                                       trait = "Rate", 
#                                                       gen = "Both", 
#                                                       fruit1 = "Strawberry", 
#                                                       fruit2 = "Cherry"))
# 
# (PAIR_STRW_BLACK_RATE <- plot_PairwisePOP_residuals(dataset = data_PERF_Rate, 
#                                                     trait = "Rate", 
#                                                     gen = "Both",
#                                                     fruit1 = "Blackberry", 
#                                                     fruit2 = "Strawberry"))
# 
# ################# PREFERENCE
# 
# # ###PLOT 
# (PAIR_BLACK_CHERRY_G2_PREF <- plot_PairwisePOP_residuals(dataset = data_PREF_three,
#                                   trait = "Nb_eggs", gen = "G2",
#                                   fruit1 = "Cherry", fruit2 = "Blackberry"))
# # 
# # (PAIR_CHERRY_STRAW_G2_PREF  <- plot_PairwisePOP_residuals(dataset = data_PREF_three, 
# #                                                             trait = "Nb_eggs", gen = "G2",
# #                                                             fruit1 = "Strawberry", fruit2 = "Cherry"))
# # 
# # (PAIR_STRW_BLACK_G2_PREF  <- plot_PairwisePOP_residuals(dataset = data_PREF_three, 
# #                                                           trait = "Nb_eggs", gen = "G2",
# #                                                           fruit1 = "Blackberry", fruit2 = "Strawberry"))
# 
# 
# (PAIR_BLACK_CHERRY_PREF <- plot_PairwisePOP_residuals(dataset = data_PREF_three,
#                                                       trait = "Nb_eggs", 
#                                                       gen = "Both", 
#                                                       fruit1 = "Cherry", 
#                                                       fruit2 = "Blackberry"))
# (PAIR_CHERRY_STRAW_PREF  <- plot_PairwisePOP_residuals(dataset = data_PREF_three, 
#                                                        trait = "Nb_eggs", 
#                                                        gen = "Both",
#                                                        fruit1 = "Strawberry", 
#                                                        fruit2 = "Cherry"))
# (PAIR_STRW_BLACK_PREF  <- plot_PairwisePOP_residuals(dataset = data_PREF_three, 
#                                                      trait = "Nb_eggs", 
#                                                      gen = "Both",
#                                                      fruit1 = "Blackberry", 
#                                                      fruit2 = "Strawberry"))







# legend <- lemon::g_legend(PAIR_BLACK_CHERRY_G2_PREF)
# 
# # ############## THIRD GENERATION 
# # HETEROGENEITY_PLOT_THIRD <- cowplot::ggdraw() +
# #   cowplot::draw_plot(PAIR_BLACK_CHERRY_G2_EGGS+theme(legend.position = "none", 
# #                                         plot.title = element_blank()),
# #             x = 0, y = 0.75, width = 0.25, height = 0.25) +
# #   cowplot::draw_plot(PAIR_CHERRY_STRAW_G2_EGGS+theme(legend.position = "none", 
# #                                         plot.title = element_blank()), 
# #             x = 0.30, y = 0.75, width = 0.25, height = 0.25) +
# #   cowplot::draw_plot(PAIR_STRW_BLACK_G2_EGGS+theme(legend.position = "none", 
# #                                         plot.title = element_blank()),
# #             x = 0.60, y = 0.75, width = 0.25, height = 0.25) +    
# #   #cowplot::draw_plot(legend, x = 0.85, y = 0.8, width = 0.0001, height = 0.0001) +
# #   cowplot::draw_plot(PAIR_BLACK_CHERRY_G2_RATE+theme(legend.position = "none", 
# #                                         plot.title = element_blank()),
# #             x = 0, y = 0.5, width = 0.25, height = 0.25) +
# #   cowplot::draw_plot(PAIR_CHERRY_STRAW_G2_RATE+theme(legend.position = "none", 
# #                                         plot.title = element_blank()), 
# #             x = 0.30, y = 0.5, width = 0.25, height = 0.25) +
# #   cowplot::draw_plot(PAIR_STRW_BLACK_G2_RATE+theme(legend.position = "none", 
# #                                         plot.title = element_blank()),
# #             x = 0.60, y = 0.5, width = 0.25, height = 0.25) +    
# #   cowplot::draw_plot(PAIR_BLACK_CHERRY_G2_ADULTS+theme(legend.position = "none", 
# #                                         plot.title = element_blank()),
# #             x = 0, y = 0.25, width = 0.25, height = 0.25) +
# #   cowplot::draw_plot(PAIR_CHERRY_STRAW_G2_ADULTS+theme(legend.position = "none", 
# #                                         plot.title = element_blank()), 
# #             x = 0.30, y = 0.25, width = 0.25, height = 0.25) +
# #   cowplot::draw_plot(PAIR_STRW_BLACK_G2_ADULTS+theme(legend.position = "none", 
# #                                         plot.title = element_blank()),
# #             x = 0.60, y = 0.25, width = 0.25, height = 0.25) +    
# #   cowplot::draw_plot(PAIR_BLACK_CHERRY_G2_PREF+theme(legend.position = "none", 
# #                                         plot.title = element_blank()),
# #             x = 0, y = 0, width = 0.25, height = 0.25) +
# #   cowplot::draw_plot(PAIR_CHERRY_STRAW_G2_PREF+theme(legend.position = "none", 
# #                                         plot.title = element_blank()), 
# #             x = 0.30, y = 0, width = 0.25, height = 0.25) +
# #   cowplot::draw_plot(PAIR_STRW_BLACK_G2_PREF+theme(legend.position = "none", 
# #                                         plot.title = element_blank()),
# #             x = 0.60, y = 0, width = 0.25, height = 0.25) +    
# #   cowplot::draw_plot(legend, x = 0.93, y = 0.5, width = 0.0001, height = 0.0001) 
# #   HETEROGENEITY_PLOT_THIRD
# 
# 
# # cowplot::save_plot(file =here::here("figures","Heterogeneity_Across_Pop_ThirdGeneration.pdf"),
# #                   HETEROGENEITY_PLOT_THIRD, base_height = 35/cm(1),
# #                   base_width = 35/cm(1), dpi = 1200)
# # 
# 
# 
#   
# legend <- lemon::g_legend(PAIR_CHERRY_STRAW_PREF)
# 
#   
# ############## ALL GENERATION 
# HETEROGENEITY_PLOT <- cowplot::ggdraw() +
#   cowplot::draw_plot(PAIR_BLACK_CHERRY_PREF+theme(legend.position = "none", 
#                                         plot.title = element_blank()),
#             x = 0, y = 0.66, width = 0.28, height = 0.28) +
#   cowplot::draw_plot(PAIR_CHERRY_STRAW_PREF+theme(legend.position = "none", 
#                                         plot.title = element_blank()), 
#             x = 0.30, y = 0.66, width = 0.28, height = 0.28) +
#   cowplot::draw_plot(PAIR_STRW_BLACK_PREF+theme(legend.position = "none", 
#                                         plot.title = element_blank()),
#             x = 0.60, y = 0.66, width = 0.28, height = 0.28) +    
#   #cowplot::draw_plot(legend, x = 0.85, y = 0.8, width = 0.0001, height = 0.0001) +
#   cowplot::draw_plot(PAIR_BLACK_CHERRY_EGGS+theme(legend.position = "none", 
#                                         plot.title = element_blank()),
#             x = 0, y = 0.33, width = 0.28, height = 0.28) +
#   cowplot::draw_plot(PAIR_CHERRY_STRAW_EGGS+theme(legend.position = "none", 
#                                         plot.title = element_blank()), 
#             x = 0.30, y = 0.33, width = 0.28, height = 0.28) +
#   cowplot::draw_plot(PAIR_STRW_BLACK_EGGS+theme(legend.position = "none", 
#                                         plot.title = element_blank()),
#             x = 0.60, y = 0.33, width = 0.28, height = 0.28) +    
#   cowplot::draw_plot(PAIR_BLACK_CHERRY_RATE+theme(legend.position = "none", 
#                                         plot.title = element_blank()),
#             x = 0, y = 0, width = 0.28, height = 0.28) +
#   cowplot::draw_plot(PAIR_CHERRY_STRAW_RATE+theme(legend.position = "none", 
#                                         plot.title = element_blank()), 
#             x = 0.30, y = 0, width = 0.28, height = 0.28) +
#   cowplot::draw_plot(PAIR_STRW_BLACK_RATE+theme(legend.position = "none", 
#                                         plot.title = element_blank()),
#             x = 0.60, y = 0, width = 0.28, height = 0.28) +    
#   cowplot::draw_plot(legend, x = 0.93, y = 0.5, width = 0.0001, height = 0.0001)  + 
#   cowplot::draw_plot_label(c("Oviposition preference", 
#                              "Oviposition stimulation", 
#                              "Offpsring performance"),  
#                   x = c(0.48,0.48,0.48), 
#                   y = c(0.96, 0.63, 0.29), 
#                   hjust = c(0.5,0.5,0.5), 
#                   vjust = c(0.5,0.5,0.5),
#                   size = 16) 
# 
# 
# 
# 
# 
# 
# # # 
# # cowplot::save_plot(file =here::here("figures","Heterogeneity_Across_Pop.pdf"),
# #                   HETEROGENEITY_PLOT, base_height = 30/cm(1),
# #                   base_width = 30/cm(1), dpi = 1200)
# 

9 PLOT POP EFFECT - HETEROGENEITY

# library("dplyr")
PPMR_ALL_Pref_Stim_G0 <- plot_pairwise_meanresiduals_popeffect(dataset = data_PERF, 
                                        formula="log(Nb_eggs+1) ~ Test_environment " ,
                                        gen = "G0" , 
                                        formula_Blanquart="log(Nb_eggs+1) ~  Test_environment + Population + SA + Test_environment:Original_environment",
                                        xaxis_labelprint = "Oviposition stimulation in sympatry",
                                        yaxis_labelprint = "Oviposition stimulation in allopatry")
## [1] "Converting test_environment column into a factor"
## [1] "Converting original_environment column into a factor"
PPMR_ALL_Pref_Stim_G2 <- plot_pairwise_meanresiduals_popeffect(dataset = data_PERF, 
                                        formula="log(Nb_eggs+1) ~ Test_environment " ,
                                        formula_Blanquart="log(Nb_eggs+1) ~  Test_environment + Population + SA + Test_environment:Original_environment",
                                        gen = "G2" , 
                                        xaxis_labelprint = "Oviposition stimulation in sympatry",
                                        yaxis_labelprint = "Oviposition stimulation in allopatry")
## [1] "Converting test_environment column into a factor"
## [1] "Converting original_environment column into a factor"
## [1] "Populations that do not have measures in sympatry have been removed"
PPMR_ALL_Pref_G0 <- plot_pairwise_meanresiduals_popeffect(dataset = data_PREF_three, 
                                        formula="log(Nb_eggs+1) ~ Test_environment " ,
                                        formula_Blanquart="log(Nb_eggs+1) ~ Test_environment + Population + SA + 
                    Test_environment:Original_environment + BoxID",
                                        gen = "G0" , 
                                        xaxis_labelprint = "Oviposition preference in sympatry",
                                        yaxis_labelprint = "Oviposition preference in allopatry")
## [1] "Converting test_environment column into a factor"
## [1] "Converting original_environment column into a factor"
PPMR_ALL_Pref_G2 <- plot_pairwise_meanresiduals_popeffect(dataset = data_PREF_three, 
                                        formula="log(Nb_eggs+1) ~ Test_environment " ,
                                        formula_Blanquart="log(Nb_eggs+1) ~ Test_environment + Population + SA + 
                    Test_environment:Original_environment + BoxID",
                                        gen = "G2" , 
                                        xaxis_labelprint = "Oviposition preference in sympatry",
                                        yaxis_labelprint = "Oviposition preference in allopatry")
## [1] "Converting test_environment column into a factor"
## [1] "Converting original_environment column into a factor"
PPMR_ALL_Perf_G0 <- plot_pairwise_meanresiduals_popeffect(dataset = data_PERF_Rate, 
                                        formula="asin(sqrt(Rate)) ~ Test_environment + log(Nb_eggs)" ,
                                        gen = "G0" , 
                                        formula_Blanquart="asin(sqrt(Rate)) ~  Test_environment + Population + SA + log(Nb_eggs) + Test_environment:Original_environment",
                                        xaxis_labelprint = "Offspring performance in sympatry",
                                        yaxis_labelprint = "Offspring performance in allopatry")
## [1] "Converting test_environment column into a factor"
## [1] "Converting original_environment column into a factor"
PPMR_ALL_Perf_G2 <- plot_pairwise_meanresiduals_popeffect(dataset = data_PERF_Rate, 
                                        formula="asin(sqrt(Rate)) ~ Test_environment + log(Nb_eggs)" ,
                                        gen = "G2" ,
                                        formula_Blanquart="asin(sqrt(Rate)) ~  Test_environment + Population + SA + log(Nb_eggs) + Test_environment:Original_environment",
                                        xaxis_labelprint = "Offspring performance in sympatry",
                                        yaxis_labelprint = "Offspring performance in allopatry")
## [1] "Converting test_environment column into a factor"
## [1] "Converting original_environment column into a factor"
## [1] "Populations that do not have measures in sympatry have been removed"
legend <- lemon::g_legend(PPMR_ALL_Perf_G2)


PPMR_ALL <- cowplot::ggdraw() + 
  cowplot::draw_plot(PPMR_ALL_Pref_G0 + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.01, y = 0.5, width = 0.25, height = 0.45) + 
  cowplot::draw_plot( PPMR_ALL_Pref_Stim_G0+ theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.31, y = 0.5, width = 0.25, height = 0.45) + 
  cowplot::draw_plot(PPMR_ALL_Perf_G0 + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.61, y = 0.5, width = 0.25, height = 0.45) + 
  cowplot::draw_plot(legend, x = 0.93, y = 0.5, width = 0.0001, height = 0.0001) + 
  cowplot::draw_plot(PPMR_ALL_Pref_G2 + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.01, y = 0, width = 0.25, height = 0.45) + 
  cowplot::draw_plot(PPMR_ALL_Pref_Stim_G2 + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.31, y = 0, width = 0.25, height = 0.45) + 
  cowplot::draw_plot(PPMR_ALL_Perf_G2 + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.61, y = 0, width = 0.25, height = 0.45) + 
  cowplot::draw_plot_label(c("Generation G0/G1", "A", "B", "C", " ",
                    "Generation G2/G3", "D", "E", "F", " "),  
                  x = c(0.4, 0.01, 0.30, 0.61, 0.92, 0.4, 0.01, 0.3, 0.61, 0.92), 
                  y = c(1, 0.98, 0.98, 0.98, 0.98, 0.5, 0.48, 0.48, 0.48, 0.48), 
                  hjust = c(0,0,0,0,0,0,0,0,0,0), 
                  vjust = c(1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5),
                  size = 16) 
 PPMR_ALL

# 
# cowplot::save_plot(file =here::here("figures", "POP_effect_ALL.pdf"),
#                   PPMR_ALL,
#                    base_height = 20/cm(1), base_width = 36/cm(1), dpi = 610)
# 

10 RELATION BETWEEN TRAITS

PLOT_Blackberry <- plot_RelationTraits_residuals(gen = "G2", 
                                                fruit = "Blackberry", 
                                                trait2 = "Preference")
PLOT_Cherry <- plot_RelationTraits_residuals(gen = "G2", 
                                                fruit = "Cherry", 
                                                trait2 = "Preference")
PLOT_Strawberry <- plot_RelationTraits_residuals(gen = "G2", 
                                                fruit = "Strawberry", 
                                                trait2 = "Preference")

PLOT_Blackberry_G0 <- plot_RelationTraits_residuals(gen = "G0", 
                                                fruit = "Blackberry", 
                                                trait2 = "Preference")
PLOT_Cherry_G0 <- plot_RelationTraits_residuals(gen = "G0", 
                                                fruit = "Cherry", 
                                                trait2 = "Preference")
PLOT_Strawberry_G0 <- plot_RelationTraits_residuals(gen = "G0", 
                                                fruit = "Strawberry", 
                                                trait2 = "Preference")

legend <- lemon::g_legend(PLOT_Strawberry_G0)



## ALL GENERATIONS 
 RELATION_TRAITS_G0G2 <- cowplot::ggdraw() + 
  cowplot::draw_plot(PLOT_Blackberry_G0 + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.01, y = 0.5, width = 0.24, height = 0.43) + 
  cowplot::draw_plot(PLOT_Cherry_G0 + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.31, y = 0.5, width = 0.24, height = 0.43) + 
  cowplot::draw_plot(PLOT_Strawberry_G0 + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.61, y = 0.5, width = 0.24, height = 0.43) + 
  cowplot::draw_plot(legend, x = 0.93, y = 0.5, width = 0.0001, height = 0.0001) + 
  cowplot::draw_plot(PLOT_Blackberry + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.01, y = 0, width = 0.24, height = 0.43) + 
  cowplot::draw_plot(PLOT_Cherry + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.31, y = 0, width = 0.24, height = 0.43) + 
  cowplot::draw_plot(PLOT_Strawberry + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.61, y = 0, width = 0.24, height = 0.43) + 
  cowplot::draw_plot_label(c("Generation: G0/G1", "A", "B", "C", " ",
                    "Generation: G2/G3", "D", "E", "F", " "),  
                  x = c(0.4, 0.01, 0.30, 0.61, 0.92, 0.4, 0.01, 0.3, 0.61, 0.92), 
                  y = c(1, 0.98, 0.98, 0.98, 0.98, 0.5, 0.48, 0.48, 0.48, 0.48), 
                  hjust = c(0,0,0,0,0,0,0,0,0,0), 
                  vjust = c(1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5),
                  size = 16) 
 RELATION_TRAITS_G0G2

# 
# cowplot::save_plot(file =here::here("figures", "RELATION_TRAITS_G0G2.pdf"),
#                   RELATION_TRAITS_G0G2,
#                   base_height = 18/cm(1), base_width = 34/cm(1), dpi = 610)

#  


##### With stimulation 
 PLOT_Blackberry_Stim <- plot_RelationTraits_residuals(gen = "G2", 
                                                fruit = "Blackberry", 
                                                trait2 = "Stimulation")
PLOT_Cherry_Stim <- plot_RelationTraits_residuals(gen = "G2", 
                                                fruit = "Cherry", 
                                                trait2 = "Stimulation")
PLOT_Strawberry_Stim <- plot_RelationTraits_residuals(gen = "G2", 
                                                fruit = "Strawberry", 
                                                trait2 = "Stimulation")
PLOT_Blackberry_Stim_G0 <- plot_RelationTraits_residuals(gen = "G0", 
                                                fruit = "Blackberry", 
                                                trait2 = "Stimulation")
PLOT_Cherry_Stim_G0 <- plot_RelationTraits_residuals(gen = "G0", 
                                                fruit = "Cherry", 
                                                trait2 = "Stimulation")
PLOT_Strawberry_Stim_G0 <- plot_RelationTraits_residuals(gen = "G0", 
                                                fruit = "Strawberry", 
                                                trait2 = "Stimulation")
 
 
 ## ALL GENERATIONS 
 RELATION_TRAITS_Stimulation_G0G2 <- cowplot::ggdraw() + 
  cowplot::draw_plot(PLOT_Blackberry_Stim_G0 + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.01, y = 0.5, width = 0.24, height = 0.43) + 
  cowplot::draw_plot(PLOT_Cherry_Stim_G0 + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.31, y = 0.5, width = 0.24, height = 0.43) + 
  cowplot::draw_plot(PLOT_Strawberry_Stim_G0 + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.61, y = 0.5, width = 0.24, height = 0.43) + 
  cowplot::draw_plot(legend, x = 0.93, y = 0.5, width = 0.0001, height = 0.0001) + 
  cowplot::draw_plot(PLOT_Blackberry_Stim + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.01, y = 0, width = 0.24, height = 0.43) + 
  cowplot::draw_plot(PLOT_Cherry_Stim + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.31, y = 0, width = 0.24, height = 0.43) + 
  cowplot::draw_plot(PLOT_Strawberry_Stim + theme(legend.position = "none", 
                                        plot.title = element_blank()), 
            x = 0.61, y = 0, width = 0.24, height = 0.43) + 
  cowplot::draw_plot_label(c("Generation: G0/G1", "A", "B", "C", " ",
                    "Generation: G2/G3", "D", "E", "F", " "),  
                  x = c(0.4, 0.01, 0.30, 0.61, 0.92, 0.4, 0.01, 0.3, 0.61, 0.92), 
                  y = c(1, 0.98, 0.98, 0.98, 0.98, 0.5, 0.48, 0.48, 0.48, 0.48), 
                  hjust = c(0,0,0,0,0,0,0,0,0,0), 
                  vjust = c(1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5),
                  size = 16) 
 RELATION_TRAITS_Stimulation_G0G2

# 
# cowplot::save_plot(file =here::here("figures", "RELATION_TRAITS_Stimulation_G0G2.pdf"),
#                   RELATION_TRAITS_Stimulation_G0G2,
#                   base_height = 18/cm(1), base_width = 34/cm(1), dpi = 610)